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Abstract 

We use high order finite difference methods to solve the wave equation 
in the second order form. The spatial discretization is performed by finite 
difference operators satisfying a summation-by-parts property. The focus 
of this work is on the numerical treatment of non-conforming grid inter¬ 
faces. The interface conditions are imposed weakly by the simultaneous 
approximation term technique in combination with interface operators, 
which move the discrete solutions between the grids on the interface. In 
particular, we consider interpolation operators and projection operators. 

A norm-compatibility condition, which leads to stability for first order 
hyperbolic systems, does not suffice for second order wave equations. An 
extra constraint on the interface operators must be satisfied to derive an 
energy estimate for stability. We carry out eigenvalue analyses to inves¬ 
tigate the additional constraint and how it is related to stability, and 
find that the projection operators have better stability properties than 
the interpolation operators. In addition, a truncation error analysis is 
performed to study the convergence property of the numerical schemes. 

In the numerical experiments, the stability and accuracy properties of 
the numerical schemes are further explored, and the practical usefulness 
of non-conforming grid interfaces is presented and discussed in two effi¬ 
ciency studies. 

Keywords: Second order wave equation. Finite difference method, SBP- 
SAT, Non-conforming grid interface, Interpolation, Coupling 
AMS subject classifications: 65M06, 65M12, 65MI5 


1 Introduction 

For wave propagation problems, the computational domain is often large com¬ 
pared with the wavelength, and waves travel for a long time. It has been shown 
that high order accurate discretizations in time and space are more efficient to 
solve these problems on smooth domains [21 [U. Although it is straightforward 
to derive high order finite difference schemes in the interior of the computational 
domain, it is challenging to derive the boundary closures in a stable and accu¬ 
rate way. For long time simulations, it is also desirable that the discretization 
is strictly stable [HI pp. 129]. A successful candidate of high order finite dif¬ 
ferences is the summation-by-parts simultaneously approximation term (SBP- 
SAT) method [3l|26]. An SBP operator [13] approximates a spatial derivative, 
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and mimics integration-by-parts via its associated norm. The SAT method [5] 
is used to impose boundary conditions and grid interface conditions weakly. 

Traditionally, the wave equation is written as a first order hyperbolic system, 
and is then solved by the well-developed methods for such systems. However, 
there are various drawbacks in doing so [12]. Therefore, it is desirable to solve 
directly the wave equation in the second order form. In [^, an SBP-SAT 
method for the wave equation in the second order form is developed and the 
numerical treatment of conforming grid interfaces is discussed. Stability of the 
numerical scheme is proved by the energy method, and the convergence property 
is investigated in the numerical experiments. 

For a wave that travels in an inhomogeneous medium, the wave speed varies 
in space. Since the wavelength is proportional to the wave speed, a reduction 
in the wave speed confined to a subset of the physical domain yields a wave 
with a shorter wavelength localized in that subset. In [H], the accuracy of a 
numerical solution to a Cauchy problem is stated in terms of the number of grid 
points per wavelength. For computational efficiency it is important that a fine 
mesh is used in the subset that constitutes the slower media, and a coarse mesh 
elsewhere. 

To achieve this, one approach is to partition the computational domain into 
blocks, where the mesh sizes are constant in each block but differ in differ¬ 
ent blocks. In more than one space dimension, the partition results in non- 
conforming grid interfaces with hanging nodes. Suitable interface conditions 
are then imposed to couple adjacent mesh blocks and yield a well-posed prob¬ 
lem. Many techniques for the numerical treatment of interface conditions have 
been proposed. In m, an energy conserving discretization of the elastic wave 
equation in the second order formulation is presented. The finite difference op¬ 
erators satisfy the principle of SBP, and the grid interface with a 1:2 refinement 
ratio is handled by ghost points. Stability is proved by the energy method, but 
the convergence rate of the numerical scheme is limited to two. 

In [16] the norm-compatible interpolation operators are constructed to han¬ 
dle non-conforming grid interfaces, and the Euler equations are used as the 
model problem. The norm-compatibility condition leads to an energy estimate 
for first order hyperbolic systems and the Schrddinger equation [22] . The inter¬ 
polation error is of the same magnitude as the error due to the derivative ap¬ 
proximations by the SBP operators. In the numerical experiments, it is shown 
that the convergence rate of the numerical scheme applied to the Euler equa¬ 
tions is not lowered by using the interpolation operators, compared with the 
case with only conforming grid interfaces. To use the interpolation operators 
presented in [16] , the mesh refinement ratio is fixed to 1 : 2 and the mesh blocks 
must be conforming. 

Recently, a general purpose methodology for coupling mesh blocks with non- 
conforming interfaces was developed in m- This technique uses projection 
operators to move a discrete finite difference solution to piecewise polynomial 
functions in a subspace of a Hilbert space where the coupling is done. The 
wave equation in the first order system formulation is used as the model prob¬ 
lem. Stability is proved by the energy method, and the numerical experiments 
demonstrate that the convergence rate is the same as if conforming grid inter¬ 
faces were used. The projection operators allow for a very flexible configuration 
of meshing in the sense that the interfaces as well as the mesh blocks do not 
need to be conforming. Similarly to the interpolation operators, the projec- 
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tion operators satisfy norm-compatibility conditions that are essential for the 
stability proof to hold. 

In this paper, we focus on the numerical treatment of non-conforming grid 
interfaces for the wave equation in the second order form in the framework of the 
SBP-SAT methodology. In particular, the stability and accuracy properties are 
investigated. We have found that in contrast to first order hyperbolic systems 
for which the norm-compatibility condition leads to stability, an extra condition 
on the interface operator is needed to derive an energy estimate for the second 
order wave equation. This condition is satisfied for the second and fourth order 
accurate interpolation operators constructed in |I4) and the projection operators 
in [To]. We prove stability by the energy method for those cases. For higher 
order accurate schemes the extra condition is not satisfied, and we cannot prove 
stability. With an eigenvalue analysis, we have found that the violation of the 
stability condition is very weak for the projection operators, and in all the 
numerical experiments we have conducted no unphysical growth is observed 
for the schemes with the projection operators. In certain cases the SBP-SAT 
schemes with the sixth and eighth order accurate interpolation operators are 
not stable. 

Local mesh refinement reduces the number of grid points significantly in 
computations. To achieve full efficiency, the numerical scheme must also be 
accurate enough. It is desirable that the convergence rate is not depressed by 
using non-conforming grid interfaces. Even though this is in most cases true 
for first order hyperbolic systems, the situation for second order equations is 
less favourable. By a truncation error analysis, we show that the truncation 
error near the edge of the non-conforming grid interfaces is two orders larger 
than that with conforming grid interfaces. The large truncation error is only 
localized at a few grid points in a two dimensional space, and its effect to the 
convergence rate may be weakened. In fact, the numerical experiments show 
that the convergence rate with a non-conforming grid interface is only one order 
lower than the corresponding case with a conforming grid interface. In addition, 
an efficiency study is carried out by a comparison of the numerical schemes with 
interpolation operators and projection operators. We have found that in certain 
cases it is beneficial to use non-conforming grid interfaces, albeit the accuracy 
reduction. 

The structure of this paper is as follows. In ^ the SBP-SAT methodology is 
introduced. We then discuss stability and accuracy properties of the numerical 
coupling based on interpolation operators in ^ and projection operators in 0 
Numerical experiments are carried out in ^ including the eigenvalue analyses 
for stability, convergence verifications for accuracy and studies on computational 
efficiency. We conclude and mention future work in ^ 


2 Preliminaries 

We begin with the preliminaries that will be used in the discussion of the SBP- 
SAT method. Let wi{x) and W 2 {x) be two real-valued functions in L2[0,1]. The 
inner product is defined by {wi,W 2 ) = Jq wiW 2 dx. The corresponding norm is 
llwilP = (wi,ryi). The computational domain [0,1] is discretized by iV + I 
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equidistant grid points 

Xi = ih, i = 0,1, ■ ■ ■ , N, where ^ ~ 

With any fixed N, a grid function can be represented by a vector and an operator 
can be represented by a matrix. Throughout this paper, we use an operator and 
a matrix interchangeably when there is no ambiguity. 

2.1 The SBP operators 

We need the following definitions: 

Definition 2.1 A difference operator Di = H~^Q approximating first deriva¬ 
tive d/dx is a narrow diagonal first derivative SBP operator, if H is diagonal 
and positive definite, Q + = B = diag{—l, 0,..., 0,1) and the interior stencil 

width of Di is minimal. 

Definition 2.2 A difference operator -\-B^^'iS) approximat¬ 

ing second derivative d/dx{b(x)d/dx) with b(x) > 0 is a narrow diagonal vari¬ 
able coefficient second derivative SBP operator, if PI is diagonal and positive def¬ 
inite, is symmetric positive semidefinite, B^^^ = diag{—b{xo), 0 ,..., 0 , b{xN)), 
the first and last row of S is a one sided approximation of djdx at the boundary 
and the interior stencil width of is minimal. 

The diagonal positive definite matrix PI defines the SBP norm, and it has 
the interior weight h and special boundary weights. To solve the wave equation 
on a curvilinear grid, even if the original equation has only constant coefficients, 
the transformed equation has second derivative terms with variable coefficients, 
and mixed derivative terms. Therefore, it is important that Di and are 
based on the same norm H. In addition, Di and must be compatible for the 
energy method to be applicable for proving stability. Compatibility m means 
that can be written as where is symmetric 

positive semidefinite. In this case, is essentially equal to applying Di twice 
plus a small dissipative term. 

The accuracy of the SBP operators are often termed as 2p, meaning that 
the approximation error of Di and is Ofh'^^) in the interior, while near 
the boundary the approximation error increases to 0{hP). For S' ~ ^ at the 
boundary, the approximation error is 0{hP~^^). 2p*^ order accurate SBP oper¬ 
ators Di are constructed in for p = 1,2, 3,4, and in [T5] for p = 5. In [H] , 
2p*^ order accurate are constructed for p = 1,2, 3, and they are compatible 
with Di constructed in The construction of requires to solve a large 
system of nonlinear equations, which makes it very involved when the accuracy 
order is high. In certain cases, we only need an SBP operator D 2 that approxi¬ 
mates second derivative d'^ jdx^, and these operators are constructed in [12] for 
p= 1, 2, 3,4 and in [T5j for p = 5. 

The following lemma, which is often referred to as the borrowing trick m, 
describes an important property of the SBP operator D 2 constructed in [Ulin]: 
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2p 

2 

4 

6 

8 

10 

Oi2p 

0.4 

0.2508560249 

0.1878715026 

0.0015782259 

0.0351202265 


Table 1: a 2 p values 


Lemma 2.3 The matrix M in D 2 constructed in m can be written as 

M = ha2p{BSf BS + M, 

where M is symmetric positive semidefinite and a 2 p is a constant independent 
of h. The values of a 2 p are listed in Table{J\for different accuracy orders. 

A similar lemma for constructed in [T3] is presented in m- 
There are other types of SBP operators as well. It is possible to increase 
the accuracy by using a block norm H, meaning that H is diagonal in the 
interior and has a symmetric-block-structured boundary closure. In this case, 
the approximation error in the interior and near the boundary are and 

0{h^P~^), respectively. The drawback of the block norm SBP operators is that 
the energy method is not applicable to prove stability for the wave equation with 
variable coefficients and/or solved on a curvilinear grid, and unphysical growth 
is indeed observed in numerical experiments. This has limited the practical 
usage of the block norm SBP operators. In m, artificial dissipation is added 
to stablize the scheme. 


3 Non—conforming grid interfaces handled by in¬ 

terpolation operators 

The numerical coupling of conforming grid interfaces by the SAT method for 
the wave equation is discussed in m- Our aim in this section is to generalize 
the scheme to also couple non-conforming grid interfaces by using the interpo¬ 
lation operators constructed in |16] . To begin with, we consider the two-block- 
structured mesh 0 shown in Figure [Tal a coarse mesh 0/, in the left block with 
nxL X nyL grid points and a fine mesh fin in the right block with UxR x Uyn 
grid points. The equality UyR = 2nyr — 1 yields a 1:2 mesh refinement ratio 
across the interface. 

In [16], the interpolation operators are denoted by If 2 C and Ic 2 F, where the 
subscripts F2C and C2F refer to fine to coarse and coarse to fine, respectively. 
Though only interpolation operators for the grid interface with mesh refinement 
ratio 1:2 are reported, the construction of the interpolation operators with some 
other refinement ratios (1:4, 1:8, •••) is possible by the same technique. To 
handle a multi-block-structured mesh shown in Figure M even though the 
interface between the left-up block and the left-down block is conforming, it is 
necessary to treat it as an interface to make the energy method applicable to 
prove stability. In other words, the mesh blocks must be conforming. 

It is important that the interpolation operators preserve the SBP property. 
To this end, the following norm-compatibility condition is essential: 

HyFlIc2F = {HyLlF2c)"'", 


( 1 ) 
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(a) Two mesh blocks (b) Four mesh blocks 

Figure 1: Non-conforming grid interfaces 


where Hy^ and Hyfi are SBP norms in the left block and right block, both in 
the y direction. 

The interpolation operators do not interpolate exactly, instead they mimic 
the accuracy properties of the diagonal norm SBP operators. The interpolation 
error is 0{h‘^P) in the interior of the interface and 0{hP) near the edge of the 
interface. We call them order accurate interpolation operators, and when 

used together with the order accurate SBP operators the scheme is also 

termed as 2p*^ order accurate, though the truncation error of the semidiscretized 
equation may not be 0{h?^) or 0{hP). Here, h is used to denote the magnitude 
of the mesh sizes for the sake of a simplified notation, though at most four 
different mesh sizes could be present in the mesh il. 

With the above accuracy requirement and the norm-compatibility condition 
o, 2 pth Qi-(ter accurate interpolation operators are constructed for p = 1,2,3 
and 4 in [16]. Though the accuracy is reduced near the edge of the interface, the 
number of grid points with the large interpolation error 0{hP) is independent 
of h. 

As will be seen later, when the interpolation operators are used to solve 
the wave equation an extra condition is posed on the interpolation operators in 
order to apply the energy method to prove stability: 

:= HyL^IyL — If2cIc2f) > 0 , Er := HyR^IyR — Ic2fIf2c) > 0 , ( 2 ) 

where S > 0 means that the matrix S is symmetric positive semidefinite, and 
lyL, lyR StTe identity matrices. It is straightforward to show that and 'E.r 
are symmetric, but the positive semidefiniteness is not a built-in constraint 
in the construction process of the interpolation operators in m- In ^ an 
eigenvalue analysis is performed to show that is satisfied for the second and 
fourth order accurate interpolation operators. Negative eigenvalues are present 
with the sixth and eighth order accurate interpolation operators. However, the 
scheme with sixth other accurate interpolation operators seems stable in some 
of the numerical experiments conducted in ^ indicating that ([5]) is a sufficient 
but not necessary condition for stability. 
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Matrix 

Size 

Nonzero 

Eqr 

“^xL ^ ^xL 

{'^xL: '^xl) 

Eqr 

'^xR ^ ’^xR 

(1,1) 

Err 

TlxL ^ '^xR 

{nxL, 1) 

Err 

'^xR ^ f^xL 

(1; '^xR) 


Table 2: Matrices 
that are used to 
pick up solutions 
on the interface. 


3.1 The wave equation with a non—conforming grid inter¬ 
face 

The wave equation in the second order form in two space dimensions is 

Utt = Uxx + Uyy, (3) 

where — oo < x < oo, 0 < y < 1 and 0 < t < tf. We assume that the 
initial conditions and boundary conditions are compatible smooth functions 
with compact support. As a consequence, the true solution U is also smooth. 
To solve the equation on the mesh shown in Figure Hal continuity of the 
solution and continuity of first normal derivative across the grid interface are 
required. 

In the numerical coupling scheme, we frequently pick up solutions along the 
interface by using the matrices defined in the first column of Table [H In each 
of those matrices, all elements are zero except one element that is equal to one. 
The sizes along with the positions of the nonzero element are listed in the second 
and third column of Tabled Note that Elr = 

Bold letters are used to denote the operators in two space dimensions, which 
are obtained from the corresponding one dimensional operators through the 
Kronecker product = Ax ® ly and Ay = Ix ® Ay, where Ix and ly are 
identity matrices. 

Next, Equation (|31) is discretized by the SBP operators on and the grid 
interface conditions are imposed weakly by the SAT method and interpolation 
operators. The semidiscretized equation corresponding to reads: 

Utt — D 2 LU SATui + SATu 2 + SATqu, (4a) 

vtt — D 2 RV SAT^i SATj,2 SATq^, (4b) 

where 

SATui ~ - H^^S^^{Eolu— {Err ® If2c)v), 

SATu2 = —T {Err ® If2c)v), 

SATqu = ”2 {EqlSxLU— {Err ® If2c)SxRv), 

and 

SAT^X = -- H,^^S'^ji{EorV— {Err ® Ic2f)u), 

SAT^2 = —T H^^{EqrV— {Err ® Ic2f)u), 

SATqt, = - H,^^{EorSxRV— {Err ® Ic2f)SxLu). 








Here, u and v are grid functions in and flu, respectively. D 2 L and D 2 R 
are SBP operators approximating second derivatives. In Q the penalty terms 
for the boundary conditions are omitted as the focus here is the numerical 
treatment of the grid interface. Both SATux and SATu -2 impose weakly the 
continuity of the solution across the grid interface. The term in SATui 
makes the semidiscretization symmetric with respect to the SBP norms. The 
penalty parameter r in SAT^-z controls the strength of the weak enforcement, 
and its value is determined by the energy method. The penalty term SATqu 
imposes weakly the continuity of the first normal derivative across the grid 
interface. The penalty terms in Equation (I4bll are constructed in a similar way. 
In the following, the stability of (HJ is proved by the energy method. 

3.1.1 Stability analysis by the energy method 

Multiplying Equation (l4all by uJHl and Equation (l4bl) by vJHr, and using 
the equality Elr = and relation ([T]), 



(5) 


—v'^ 0 Hyii)v - v'^{EorSxr 0 HyR)v - tv’^{Eor 0 HyR)v 

— {S'^rErr 0 HyRlF2c)v + v'^{S'^ rErr 0 HyRlc2F)u 
+ 2 tU^[Err 0 HyRlF2c)v). 

Next, the following equality is obtained by moving all terms on the right hand 
side of to the left 



( 6 ) 


where 


=ujHLUt + vjHn.vt 

+ U^{M^r 0 HyR)u - U^{EQRSa:L ® HyR)u + TU^ {Eqr 0 HyR)u 
+ v^{M^r 0 HyR)v + v'^{EorSxR 0 HyR)v + tv'^{Eor 0 HyR)v (7) 
+ {S^rErr 0 HyRlF2c)v — {S^rErr 0 E[yRlc2F)u 
— 2 tU^{Err 0 HyRlF2c)v. 

Then, an energy estimate exists if E^ > 0 for any u and v, and in this case 
E^ is a discrete energy of the semidiscretized equation (|3]). To achieve this, 
we need relation to obtain 

(Eql 0 HyR)u > vF{E qr 0 E[yRlF2cIC2F)u 


= {Eqr 0 Ic2FHyRlc2F)u 

= {{EqL ® Ic2f)u)^{I xL ® HyR){{EoR 0 Ic2f)u). 


It then follows for u: 
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and for v. 

{EQR®HyR)v > -jv'^{EQR®HyB)v + -{{EQR®lF2c)v)'^{IxR®HyL){{EQR®lF2c)v). 


In addition, Lemma [Q gives: 

MxL = MxL + hxLa2p{EoLSxL)^ {EqlSxl), 

Mxr = Mxr + hxRa2p{EoFSxR)^ {EqrSxr), 

where Mxl and Mxr are symmetric positive semidefinite matrices. 
Next, we plug in (l2.dL (jS]) and ([S]) to ([7]), and obtain 


(9) 

( 10 ) 


pW 


— Qi + Q 2 + Qs, 


where 

Qi =ujHLUt + vfHRVt + u^{MxL ® HyL)u + v^{Mxr ® HyR)v, 
Q 2 — ^xL^2p^EqIjSxL'^^ {^xL ^ EyF^i^EQRSxL^') 

— (EqlSxlu)^ {IxL ® HyF){u — {Elf ® If2c)v) 

+ 2 ~ {Elr ® If2c)v)^{I xL ® HyF){u — (Elf (g) If2c)v), 

Qs —hxRO'2p{EoRSxRv)^{IxR (g> HyF){EoRSxRv) 

— (EorSxRv)^ {IxR 0 HyF){{EFL 0 Ic2f)u — v) 

+ -^{{EfL 0 Ic 2 f)u — v)'^{IxR 0 HyF){{EFL 0 Ic 2 f)u — v). 


Since Mxl and MxR are positive semidefinite, we have Qi > 0. To ensure 
Q 2 > 0 and Qs > 0, we need 


2 \/ha,La 2 p\A 72 > 1 and 2 ^hxRa 2 p\/T^ > 1. 


=> T > max 


1 


2cX2phxL 


1 

2a2ph 



In the mesh shown in Figure Hal the mesh refinement ratio across the grid 
interface is 1 : 2, i.e. hxR = \hxL- Thus, 


T > 


1 

^2phxL 


( 11 ) 


is the condition for E^ to be a discrete energy. Therefore, an energy estimate 
is obtained if (HI, (I2|) and (fTTll hold. 

As shown above, in order to apply the energy method must be satisfied, 
which means that the energy estimate is only valid for the second and fourth 
order accurate schemes. The same interpolation operators are used for the 
Schrodinger equation in [35] but ([31 is not needed for stability. For the Euler 
equations, (ED is not needed when the standard SBP operators are used for the 
discretization, but needed when upwind SBP operators are used [16]. 
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3.1.2 Convergence rate 

We discuss the accuracy properties of (|3]) by analyzing the truncation error 
of (Hall . Equation llbl can be analyzed in a similar way. The approximation 
error of the SBP operator £>2i is 0{h?^) in the interior and 0{hP) near the 
interface, with the latter one being the dominant source of error. In the first 
penalty term SATux, a large interpolation error 0{hP) is located near the 
edge of the grid interface. Due to the h~^ factor in both H~^ and S^j^, the 
localized truncation error of SATui is 0{hP~^). Similarly, we find that the 
localized truncation errors of SATu 2 and SATqu are and 

respectively. Therefore, the localized truncation error of the semidiscretization 
(Ha|) is 0{hP~'^). 

In [55j, the convergence of the SBP-SAT discretization of the second order 
wave equation in one space dimension with a grid interface is analyzed. The 
result is that if the penalty parameter is chosen strictly larger than the limit 
value required for stability, the localized truncation error 0{h'P) near the grid 
interface results in an error 0{hP^^) in the solution for p = 1, and an error 
for p > 2. In other words, we gain one order in convergence if p = I 
and two orders if p > 2. 

In our case, the spatial dimension is two and there is the possibility of another 
gain in convergence. That is, the number of grid points with truncation error 
0{hP~'^) is finite and independent of h. Hence, the L2 norm of this truncation 
error is 0{hP~^), and is one order higher than the pointwise truncation error. 
Therefore, we can hope to get an extra gain in convergence comparing with the 
corresponding one dimensional case. 

By a convergence test in ^ we find that the extra gain is one order, which 
gives a total gain of two orders for p = 1 and three orders for p > 2. That is, 
the localized truncation error 0{hP~^) results in an error 0{h^) in the solution 
for p = 1, and an error for p > 2. To obtain this convergence rate, 

it is important to choose the penalty parameter strictly larger than the value 
required for stability. 

Comparing with the case of conforming grid interfaces, the convergence rate 
is one order lower. Even though a non-conforming grid interface allows for a lo¬ 
cal mesh refinement, the loss of accuracy may attenuate its efhciency in practice. 
To overcome the accuracy reduction by the non-conforming grid interfaces, we 
have tried to build interpolation operators with error 0{h?^) in the interior and 
error near the edge of the interface, based on both diagonal and block 

norm SBP operators by using the symbolic software MAPLE and the approach 
presented in [16j . However, we could not find a solution to the resulting system 
of equations. 

3.2 An extension to T—junction interfaces 

In Figure flbl the interface between the two blocks on the left is conforming. It is 
then in many cases desirable not to consider it as a grid interface, but to use the 
mesh shown in Figure [^ instead, where the interface there forms a T-junction. 
This is because SBP operators have a larger approximation error near the grid 
interface than that in the interior. The usage of redundant grid interfaces results 
in additional errors in the solution. Moreover, avoiding T-junction interfaces 
may end up in a bad partitioning of the computational domain with many un- 
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(a) 


(b) 



Figure 2: T-junction interfaces 


necessary mesh blocks. With a straightforward application of the interpolation 
operators, instability occurs around the T-junction interface intersection point. 
The reason is that near the interface intersection point in Figure [5^ the SBP 
norm in the vertical direction has the interior weights in the left domain, and 
the boundary weights in the two domains on the right. As a consequence, the 
norm-compatibility condition is violated and no energy estimate can be derived. 

In [21) . the T-junction operators are constructed to handle T-junction in¬ 
terfaces and are applied to the advection equation and the Schrodinger equation 
in the SBP-SAT framework. Stability is proved by the energy method, but it 
comes with the cost that the T-junction operators introduce an error 0{h?P) 
in the interior of the interface and 0{hP) near the edge of the interface. The 
T-junction operators can also be used together with the interpolation operators 
to handle non-conforming grid interfaces. One constraint for the T-junction 
operators is that the interface intersection point must be a grid point in all 
involved mesh blocks, for example a close-up T-junction interface in Figure [2bl 
It is not straightforward to handle the T-junction interface shown in Figure!^ 
by the same technique. 


4 Non—conforming grid interfaces handled by pro¬ 

jection operators 

In [10] , a new methodology of handling grid interfaces is introduced. In contrast 
to the interpolation operators which are based on a direct interpolation tech¬ 
nique, the new methodology is based on a projection method. The highlights 
are that there is no strict requirement on the mesh refinement ratio, and the 
mesh blocks do not need to be conforming. 

To illustrate how the projection method works, we consider again the mesh 
n shown in Figure Hal and denote in and in fl/j the grids on the 
interface. In addition, and denote the discrete finite difference solutions 
on y^ and y^. In a general setting, eight projection operators are used to move 
z^ and z^ between the two grids on the interface. Firstly, the discrete finite 
difference solution z^ is projected by a projection operator to a piecewise 
continuous function based on the grid y^. The associated mass matrix on 
y^ is diagonal positive definite since Jacobi polynomials are used as the basis 
functions. Next, the glue grid y with the mass matrix is defined as the grid 
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that consists of the grid points on both and y^. The projection operator 
Pp 2 g i® used to project the piecewise continuous function from y^ to ?/, and is 
viewed as a basis transformation between polynomial spaces. Similarly, Pp 2 f 
and Pg 2 p are projection operators in the reversed direction corresponding to 
Pf 2 p and Pp 2 g, respectively. Pp 2 p, P^f, P^g and P^ 2 p are the corresponding 
projection operators for the grid y^. 

Similar to the interpolation operators, there are norm-compatibility condi¬ 
tions for the projection operators: 

HyLP^2f = {M^Phpf, HynPp% = {M^PP2pf, ( 12 ) 

and 

M^Pp% = {M^Py%f, Map^2f = {M^P^2pf- (13) 

Define 

~ P^fP^pPp 2 gPf 2 p and I^ 2 L ~ Pp 2 f Pg 2 pP^gPf 2 pi (14) 

the operator I^ 2 R moves from y^ to and I^ 2 L moves from y^ to y^. 

For the projection operators Pfl^ and Pp 2 ^ which move the discrete finite 
difference solution to the subspace of piecewise continuous functions and back, 
the projection error is 0{h?^) in the interior and 0{hP) near the edge, where 
p = 1, 2,3,4, 5. Pp 2 ^ and can be viewed as basis transformation operators 
between grids. As a consequence, the projection error of I^ 2 R ^nd I^ 2 R i® ^1®° 
0{h?P) in the interior and 0(}p) near the edge, where p = 1,2,3,4, 5. In other 
words, I^ 2 R ^^'l ^L 2 R have the same accuracy properties as the diagonal norm 
SBP operators and the interpolation operators in |16] . 

4.1 Interface treatment with the projection operators 

If no T-junction interface is present in the mesh, for example in Figure [T^ 
and [Tbl the operators I^ 2 R ^^'l ^R 2 L m (|Hll are used to impose grid interface 
conditions in the same way as the interpolation operators discussed in IJS] I^ 2 R 
and I^ 2 L satisfy an analogue of relation (IT|) up to tenth order accuracy, and 
relation m up to fourth order accuracy. The stability analysis and accuracy 
properties of the interpolation operators in ^are still valid for I^ 2 R ^^'l ^R 2 L- 

If a T-junction interface is present in the mesh, for example in Figure 
we cannot use I^ 2 R ^R 2 L m a direct way. Instability occurs around the 
junction point if on one side the SBP norm has the interior weights while on the 
other side it has the boundary weights, and no energy estimate can be obtained. 
To overcome the instability, the coupling is done on the glue grid. That is, we 
project finite difference solutions to the glue grid, compute the penalty terms 
there, and project them back to the finite difference grids. In this way, we avoid 
the instability caused by the SBP norms since the penalty terms are computed 
on the common glue grid. 

The projection technique is very flexible to handle grid interfaces in the sense 
that we are free to choose the interface structure, the mesh refinement ratio 
and the accuracy of the diagonal norm SBP operators. In m, the authors 
also couple the SBP finite difference method with the discontinuous Galerkin 
method, inspired by the relation between the discontinuous Galerkin spectral 
element method and the SBP-SAT finite difference method [S]. 
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Finally, we remark that even with the mesh in Figure [T^ where both inter¬ 
polation technique and projection technique are applicable, the interpolation 
operators Ic 2 f{If 2 c) are not the same as lL 2 Ri^R 2 L) ia (ITU) . The latter one 
has a wider stencil. In the construction procedure of these operators, one gets a 
system of linear equations after imposing stability and accuracy requirements. 
The solution of the linear system has a few free parameters. There are different 
ways to tune those free parameters. In |16| . the free parameters are used to 
minimize the coefficients of the leading interpolation error in L 2 norm, while 
in m the free parameters are used to minimize the distance between nearest 
eigenvalues of Pp 2 fPf 2 p for a finite difference grid of size 64. The choice of tun¬ 
ing free parameters has no influence on the theoretical order of accuracy, but 
may have an impact on condition ([2]) and the practical accuracy. This is studied 
in more detail in the numerical experiments in ^ 


5 Numerical experiments 

In this section, numerical experiments are performed to compare the schemes 
with the interpolation operators and the projection operators, and verify their 
stability and accuracy properties. Moreover, we also conduct two numerical 
experiments to study the efficiency of local mesh refinement by solving the wave 
equation on a domain with a complex geometry. 

The L 2 error and maximum error are computed as the norm of the difference 
between the exact solution and the numerical solution according to 

\\u^ — m®^||l2 = \J — It®®), 

||u^ — u®^||oo = max \u^ — u“|/amp, 

where d is the dimension of the equation and amp is the maximum amplitude 
of the solution. The convergence rate is computed by 

5.1 Stability study 

We begin with an eigenvalue analysis for condition (2)). The computational 
domain is a; G [—1,1] and y G [0,1] with a grid interface at a: = 0. In the left 
domain the number of grid points is 26 in both x and y directions, while in the 
right domain the number of grid points is 51 in both x and y directions. The 
mesh refinement ratio is 1 : 2 across the grid interface, which enables both the 
interpolation operators and projection operators applicable. The matrices El 
and Efi are symmetric, so they have real eigenvalues. We denote kL and ku the 
smallest eigenvalue of Sl and in ([2]), scaled by the mesh size: 

kL = m.in{eig{EL))/hyL kn = min(eig(Sfl;))//iyfl;. 

The reason for the scaling is that the elements in ElIEu are proportional to 

hyL/hyR. 

In Table 131 we list kL and kpi for the interpolation operators and the pro¬ 
jection operators in Column three and four, respectively. For the interpolation 
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-6.2- 10-1^ 
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-8.6- 10-® 
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-4.4-10-4 
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-5.0-10-3 
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-7.4-101 

-4.8- 10-4 

10 

kh 


-4.8- 10-4 


-5.0-10-3 

kR 


-1.1- 10-4 



Table 3: 10: interpolation operators, PO: projection operators. Column 3 and 
4 correspond to the numerical study of relation 0, and Column 5 and 6 corre¬ 
spond to the eigenvalue analysis of m- 


operators, @ holds for both the second and fourth order accurate cases with 
errors up to machine precision. For the second order accurate case, we can prove 
that kLjkn > 0 independent of h, because is diagonally dominant and S/j 
can be transformed to a diagonally dominant matrix without changing the signs 
of the eigenvalues. For the sixth and eighth order cases, ([2]) no longer holds. 
The difference between these two cases is that and ka are close to zero for 
the sixth order case, but far away from zero with the eighth order case. When 
increasing the number of grid points, the values of /cl and ku remain unchanged. 

For the projection operators, m also holds also for the second and fourth 
order accurate cases. For the sixth, eighth and tenth order cases, m does not 
hold anymore but the values of and k^ are close to zero, and they become 
slightly closer to zero as the mesh is refined. 

Another way of analyzing stability through numerical experiments is to write 
the semidiscretized equation 0 as a system of ordinary differential equations 

ztt ^Qz + F. (15) 

It is stable if the eigenvalues of Q are real and non-positive. Otherwise, the 
numerical scheme is unstable. We have computed the eigenvalues of Q by using 
the same mesh as for the computation of k^ and ka. All the eigenvalues are real. 
In Table [3l the largest eigenvalue of Q is shown in Column five and six for the 
schemes with the interpolation operators and projection operators, respectively. 
For the numerical scheme with the interpolation operators, it is stable for the 
second, fourth and sixth order cases, and unstable for the eighth order case. For 
the numerical scheme with the projection operators, it is stable for up to tenth 
order cases even though relation © only holds for up to fourth order scheme. 
The eigenvalue analysis indicates that m is sufficient but not necessary for 
stability. 
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5.2 Accuracy study 

In this section, the convergence of the SBP-SAT method applied to the wave 
equation m with a non-conforming grid interface is investigated. The analytical 
solution to (ISl) is manufactured, which means that a closed form is chosen and is 
used to obtain the initial and boundary data. At the outer boundaries, Dirichlet 
boundary condition is imposed weakly by the SAT method as described in [TB] . 

5.2.1 A non—conforming grid interface 

In the first numerical experiment, the computational domain is [—1,1] x [0,1] 
where a grid interface is located at x = 0,y G [0,1]. In the grid refinement 
level r = 0, the numbers of grid points in the left block and right block are 
26 X 26 and 51 x 51. The mesh sizes are halved in both x and y directions when 
r is increased by one. In this setting, the grid refinement ratio is I : 2, and 
both the interpolation operators and the projection operators are applicable for 
the numerical treatment of interface conditions. The fourth order Runge-Kutta 
method is used as the time integrator with the step size in time chosen so small 
that the temporal error is negligible compared with the spatial error. 

The manufactured solution to ([3]) is chosen to be 

U (a;, j/, t) = cos(5a; + 1) cos(5y + 2) cos(5-\/2t + 3). (16) 

The computational results are shown in Table S] where 2p and r in the 
first two columns denote the order of accuracy and the mesh refinement level, 
respectively. In Column 3, 4 and 5 the errors in L 2 norm, the convergence rates 
in L 2 norm and maximum norm are shown for the numerical schemes with the 
interpolation operators, whereas in Column 6, 7 and 8 the corresponding results 
obtained by the schemes with the projection operators are shown. 

For the scheme with the interpolation operators and time step At = O.l/i, 
the convergence rates in L 2 norm are 1, 3 and 4 for the second, fourth and 
sixth order schemes, respectively. This agrees with our accuracy discussion in 
[]3| Though stability cannot be proved by the energy method for the sixth order 
accurate case, it seems that for this particular setting the scheme is stable and 
exhibits the expected convergence rate. Instability occurs when using the eighth 
order accurate scheme. 

With the projection operators, the convergence rate in L 2 norm is one for 
the second order accurate scheme, and p + 1 for the fourth, sixth, eighth and 
tenth order accurate schemes, which agrees with the accuracy discussion in 21 
We note that though an energy estimate does not exist for the sixth, eighth and 
tenth order accurate cases, the schemes are stable and converge as expected. 
The time step is At = O.lh for 2p = 2,4, 6. With this time step, the tenth order 
accurate scheme yields slightly lower convergence rate than expected, and the 
result shown in TableUJis obtained with At = 0.05h. The eighth order accurate 
scheme is a special one, since with At = 0.05ft. it is even unstable. To obtain the 
results in Table HJ At = 0.025ft is used as the time step, which indicates that 
the eighth order accurate semidiscretized equation is stiff. Moreover, the error 
obtained with the eighth order accurate scheme is larger than the error obtained 
with the sixth order accurate scheme, except for the finest mesh refinement level. 

From Table 21 it is also observed that the errors are similar to each other 
for the schemes of the same order of accuracy with interpolation operators and 
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Table 4: Convergence of the SBP-SAT scheme for the wave equation with a grid 
interface. The interface conditions are handled by the interpolation operators 
(Column 3,4,5) and the projection operators (Column 6,7,8). 
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Table 5: Convergence for the wave equation with a T-junction interface handled 
by the projection operators. 

projection operators. 

5.2.2 A T—junction interface 

Next, we consider the computational domain [—1,1]^ that is divided into three 
mesh blocks as shown in Figure The interfaces are located at x = 0,y G 
[—1,1] and y = 0,x G [0,1]. In the grid refinement level r = 0, the numbers of 
grid points in Block 1 (left), Block 2 (right-up) and Block 3 (right-down) are 
28 X 51, 27 X 25 and 51 x 50, respectively. The mesh sizes are halved in both 
X and y directions when r is increased by one. This partitioning and meshing 
result in a highly non-conforming grid interface with a close-up shown in Figure 
[2B The interface conditions are imposed weakly by the SAT method with the 
projection operators. To test convergence, (USD is used as the analytical solution. 
The computational results are shown in table [5j 

Clearly, {p + 1)*^ convergence rate in L 2 norm is obtained for p = 2,3,4,5 
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Figure 3: An example of an acoustic time-harmonic plane wave impinging on 
a circular inclusion, where the wavelength is much smaller inside the circular 
inclusion than that outside. 


and first order convergence rate is obtained for p = 1. Here, we observe again 
that the schemes higher than fourth order accuracy are stable though no energy 
estimate can be obtained. 

5.3 Efficiency study 

In many applications, the frequencies of the present waves are given by initial 
and boundary data, and internal forcing functions. The wavelength of a wave 
is determined by the ratio between the wave speed of the material in which the 
wave is traveling and the frequency of the wave. The accuracy of a numerical 
solution can be stated in terms of how many grid points per wavelength are used 
to resolve the present waves A reduction in wave speed confined to a subset 
of the physical domain yields waves of a shorter wavelength localized in that 
subset. For accuracy it is therefore necessary to refine the grid to compensate 
for the shorter present wavelengths. For computational efficiency it is important 
that this refinement is done only in the subset that constitutes the slower media, 
since it is only in the slower media that wavelengths are reduced. As an example 
Figure [3] shows the scattering and diffraction of an acoustic time-harmonic plane 
wave impinging on a circular region of a slower material, the wavelength is seen 
to be reduced inside the circular region. 

Geometrical features of the physical domain such as a complex boundary or 
an internal cavity introduce a local radius of curvature. A small local radius of 
curvature compared with the present wavelengths can imply difficulties when 
generating a computational grid. As an example Figure Hal shows the scattering 
of an acoustic time-harmonic plane wave impinging on a circular cavity of a 
radius of curvature smaller than the wave length of the incoming and scattered 
waves. A part of the grid used to represent the wave field is shown as an inset, 
where it is seen in Figure |4b] that the quality of the grid is impaired as the grid 
spacing gets unnecessarily small close to the cavity. 

In the preceding experiments it has been verified that using interpolation and 
projection operators to patch together the computational grid in a multi-block 
fashion yields a stable discretization, the convergence rate, however, was seen to 
be reduced. In the following experiments we investigate the practical benefit of 
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(a) A circular cavity (b) A close-up of the mesh near the cavity 

Figure 4: An examples of an acoustic time-harmonic plane wave impinging on 
a circular inclusion cavity with a grid close to the cavity where the grid size in 
the azimuth direction is much smaller than that in the radial direction. 


using interpolation operators and projection operators in a region with a slower 
wave speed or a geometrical feature of a small radius of curvature albeit the order 
reduction. In particular, we will consider experiments involving acoustic waves 
impinging on a circular cavity and a circular inclusion of a differing material. 

The numerical method used to solve the acoustic wave equation in the fol¬ 
lowing experiments is based on the SBP-SAT scheme described in [^. The 
geometrical features are handled by using a multi-block strategy to decompose 
the physical domain into blocks, where each block allows for a mapping to curvi¬ 
linear coordinates. In m the blocks that constitute the domain are discretized 
by using conforming grids and patched together by the SAT method. In this pa¬ 
per we allow for non-conforming grids by implementing interpolation operators 
as well as projection operators into the handling of the multi-block interfaces. 

The following numerical experiments use two different two dimensional do¬ 
mains: 

• 2?i: An acoustic plane with a circular cavity of radius a. 

• '02'. An acoustic plane with a circular inclusion of radius a. 

The geometries of the domains are handled by decomposing each domain in a 
multi-block fashion. The blocks are then patched together to compound the 
entire domain. A detailed description of how these two grids are constructed is 
presented in the Appendix. 

5.3.1 A circular cavity 

In this numerical experiment, we consider a domain of an infinite homogeneous 
medium with a circular cavity of radius o = 1. Let a plane harmonic wave 
Ui = propagate in that domain and impinge on the cavity. A scat¬ 
tered wave Us is generated when the incident wave hits the cavity, and the 

total displacement Ui + Ug satisfies the wave equation. Homogeneous Neumann 
boundary condition is imposed at the cavity boundary. A detailed derivation of 
an analytical solution is found in |6l §7]. 
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(a) N-partitioning (b) T-partitioning 

Figure 5: Computational domain of the circular cavity experiment 


We take uj = 2 tt and c = 2.5, which give a wavelength 2.5. The compu¬ 
tational domain T>i is chosen to be the rectangular [—25.5,11.7] x [—11.7,11.7] 
and the cavity is centred at the origin. Two ways of partitioning the domain are 
considered, namely the N-partitioning and the T-partitioning shown in Figure 
l^andlKbl respectively. In the N-partitioning approach, we only use conforming 
gird interfaces and conforming mesh blocks. The cavity is surrounded by four 
blocks that constitute the square [—11.7,11.7]^, which is attached by a rectan¬ 
gular domain to the left. The numbers of grid points on each edge are shown in 
the hgure, and are chosen so that approximately 20 grid points per wavelength 
are used in the discretization. In this setting, the mesh is of bad quality since 
the mesh size near the cavity is significantly smaller than that near the outer 
boundaries. To overcome this drawback, we propose the T-partitioning where 
the cavity is surrounded by a small square block [—1.3,1.3]^. Here, all the grid 
interfaces are also conforming but a T-junction interface is present at x = —11.7 
with the intersection points marked by the dots. Again we choose the mesh size 
so that there are approximately 20 grid points per wavelength, and here it is 
only over-resolved in the small block [—1.3,1.3]^. The T-partitioning results in 
a mesh of 54903 grid points. The number of grid points with the N-partitioning 
is about doubled to 109867. 

We employ the fourth and sixth order SBP-SAT method to propagate the 
wave for ten periods, and show the recorded maximum errors in Figure and 
I6bl In both cases, the maximum error with the T-partitioning is about three 
times larger than that with the N-partitioning. The is not surprising because 
the mesh with the T-partitioning has less grid points than the mesh with the 
N-partitioning, and the corresponding scheme with the T-partitioning has one 
order lower convergence rate than that of the N-partitioning. It does not seem 
to improve the efficiency by using T-junction interfaces for this case. 

Although using T-junction interfaces introduces a larger error in the solu¬ 
tion, it could be beneficial for a problem with a more complex geometry. For 
example, if there are several cavities in the domain, an N-partitioning that only 
allows conforming mesh blocks would produce a large number of small mesh 
blocks. With a higher order accurate scheme, the boundary stencil gets wider 
and the number of grid points in each direction must be large enough in every 
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(a) Fourth order accurate scheme (b) Sixth order accurate scheme 


Figure 6 : Maximum error of the numerical experiment with the circular cavity 


mesh block. It is therefore over-resolved in those small mesh blocks and results 
in a suboptimal performance of the numerical scheme, and T-junction interfaces 
could be desirable. 


5.3.2 A circular inclusion 

Consider a circular domain of radius o = 1 embedded in an infinite surrounding 
medium of differing material with wave speed c. Let the wave speed c' of the 
circular domain be such that c' < c and let an incoming time-harmonic plane 
wave ui{x, y, t) = , 7 = ^ travel in the x-direction and impinge on the 

circular inclusion. The resulting field consists of the incoming wave u/, as well 
as the scattered and diffracted waves us and ud, respectively. The conditions 
at the interface of the circular inclusion are 


Ui +US = ud, 

g g on x'^ 


y=i, 


(17) 


where ^ denotes the normal derivative on the interface. Since d < c, the short 
wavelength occurs inside the circular domain. An analytical expression for the 
solution is given in [TJ pp. 667]. 

In the numerical experiments we take w = 27r, c = 1 and d = 1/10, which 
give a wavelength of 1 and 1/10 outside and inside the circular inclusion, re¬ 
spectively. To resolve the geometric features, the computational domain is de¬ 
composed into 10 conforming blocks as shown in Figure [T] We take the side 
length 2D — 2.6 for the square block outside the circular inclusion, and the side 
length 2d = 0.7-\/2 for the square block inside the circular inclusion. Both square 
blocks are centered at the origin. The Cartesian block [—5.9, —1.3] x [—1.3,1.3] 
is then attached to the left of this representation. Firstly, we only use conform¬ 
ing grid interfaces with the numbers of grid points in each block given in Table 
O The resolution outside the circular inclusion is about 16 and 51 points per 
wavelength in the horizontal and vertical direction, respectively. Inside the cir¬ 
cular inclusion the waves are resolved by about 10 grid points per wavelength in 
both directions. Hence, the waves are significantly over-resolved in the vertical 
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Figure 7: Computational domain of the circular inclusion experiment 
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Table 6: Number of grid points with 
conforming grid interfaces 
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Table 7: Number of grid points with 
non-conforming grid interfaces 


direction outside the inclusion, which leads to a suboptimal efficiency of the 
numerical scheme. 

To amend the over-resolution, we partition the computational domain in 
the same way as above but use non-conforming interfaces denoted by the small 
red circles in Figure [71 The non-conforming grid interfaces are handled by the 
interpolation and projection operators. The numbers of grid points in each 
block are chosen as in Table [71 Now the resolution is reduced to about 26 grid 
points per wavelength in the vertical direction outside the circular inclusion. 
The interface conditions (HZl) are imposed numerically with the SAT technique 
and at outer boundaries the exact solution is injected at all times. In [3], the 
SBP finite difference method applied to the wave equation with the injection 
method to impose the Dirichlet boundary condition is proved to be stable. 

The solution is propagated numerically for 10 temporal periods by the SBP- 
SAT method and the relative maximum error is recorded at each time step. In 
Figures |8a] and jSb] we plot the recorded relative maximum error as functions of 
time. Here we see that the errors are similar in both cases. The grid with con¬ 
forming interfaces has 46460 grid points, whereas the grid with non-conforming 
interfaces has 38710 grid points. The smallest grid size is determined by the res¬ 
olution inside the circular inclusion, for this reason the time step At = 4 x 10“'^ 
for the sixth order SBP-SAT method and At = 5 x 10“"* for the fourth order 
SBP-SAT method are the same for both grids. We conclude that even though 
the formal order of accuracy is lowered by using blocks with non-conforming 
















23 



(a) Fourth order accurate scheme 


. x10 


-ERROR WITHOUT INTERPOLATION 
-ERROR WITH PROJECTION 



4 6 

TEMPORAL PERIODS 

(b) Sixth order accurate scheme 


Figure 8: Maximum error of the numerical experiment with the circular inclu¬ 
sion 


interfaces it can be a beneficial strategy within the SBP-SAT framework when 
the physical domain contains regions that require a higher density of grid points. 
We also note that for more complex multi-block domains consisting of a larger 
number of blocks the benefits of using blocks with non-conforming grid inter¬ 
faces are expected to increase. 

In the experiment with the sixth order SBP-SAT scheme with the interpola¬ 
tion operators, the numerical solution blows up quickly, which indicates that the 
scheme is unstable. The corresponding scheme with the projection operators is 
stable. According to the stability analysis in >15.11 for the sixth order accurate 
scheme condition ([5]) holds for neither the interpolation operator nor the pro¬ 
jection operator. If the smallest eigenvalue of is non-negative, then an 

energy estimate exists that ensures stability. The smallest eigenvalue of ’^l/r 
scaled by the mesh size is in the magnitude of —10“^ for the sixth order accu¬ 
rate interpolation operators, and —10“^ for the sixth order accurate projection 
operators. The violation of (ED is much weaker for the projection operator than 
for the interpolation operator, which explains the observation in the numerical 
experiments. 


6 Conclusion and outlook 

In this work, high order accurate SBP finite difference operators are used to 
discretize the wave equation in the second order form on a block-structured 
mesh. Adjacent mesh blocks are patched together by imposing suitable interface 
conditions via the SAT technique. The emphasis is placed on the numerical 
treatment of non-conforming grid interfaces by the interpolation operators and 
projection operators, which are also compared in terms of the stability and 
accuracy properties. In contrast to first order hyperbolic equation, stability of 
the scheme for the second order wave equation introduces an extra condition on 
the numerical treatment of non-conforming grid interfaces. This condition is 
satisfied for the second and fourth order accurate cases, and an energy estimate 
is derived to ensure stability. For higher order accurate schemes, the extra 
stability condition is violated. We show by the eigenvalue analysis that the 
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violation is stronger with interpolation operators than with projection operators. 
Unphysical growths are observed in the numerical experiments with high order 
interpolation operators, whereas with projection operators the scheme is stable 
in all the experiments we have conducted. 

We have also performed a truncation error analysis and an investigation of 
the convergence property for the scheme, which indicates that the convergence 
rate is one order lower than that in the corresponding case with conforming 
grid interfaces. This is verified in numerical experiments. The efficiency studies 
show that for a problem with a very complex geometry, it could be beneficial 
to use non-conforming grid interfaces. 

For high order accurate interpolation operators and projection operators, 
there are free parameters left in the construction process. It is desirable to tune 
the free parameters so that the extra stability condition is satisfied. However, 
the resulting nonlinear problem seems difficult to solve. To overcome the accu¬ 
racy reduction, more accurate interpolation operators or new ways of imposing 
interface conditions are needed. 


Appendix 

We give a detailed description for the construction of the two grids used in 
the efficiency studies in 95.31 A general block Bi of a decomposition has four 
boundaries defined by the parametrized curves 

CiS = (XisiOtUisiO) > CiN = {XiNiO^ ViNiO) ) 

Ciw = {Xiwiv), ViWiv)), CiE = (xiE{v),yiE{v)), 

where 0<^<1,0<77<1. Cis and Cin describe one pair of opposing sides 
and Ciw and Cie the other pair. Let Visw denote the point of intersection 
between the curves Cis and Ciw e.t.c. A bijection (x, y) = Ti(^,ii) from the 
unit square S = [0,1]^ to the block Bi of the decomposition is given by the 
transfinite interpolation as 


y) = (1 — Tf)CiS + yCiN + (1 — + ^CiE 

- ^rfPiNE - C(1 - V)'PiSE - y(l - O'PiNW “ (1 “ 0(1 “ 'n)'PiSW ■ 

The unit square S is discretized by the points 

0 = jhy h = l/(^*« - 1), J = 0,..., - 1, 

Vk — — 1 / (^Nirj 1 ), k — 0 , . . . , Nirj 1 , 

where Ni^ and are integers determining the number of grid points in the 
spatial directions of the discretization of the block Bi. The corresponding grid 
points are computed as 

{xj ,yk) — Tii^^j, Tlk\ 

We now give details on how the grids in T>i —T >2 are constructed. Vi represents a 
circular cavity in an infinite surrounding media. We construct the computational 
grid such that the cavity of radius a is centered inside a square of side 2D. This 
is done by introducing four blocks B[^^-B^\ The bounding curves of block b[^^ 
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are given by 

= a - 1/^2, , 

cf^ = {D,2D^-D), 

^iw — (^~v{D — a/ \/2) — a/V2, r]{D — a/ \/2) + a/V2^ , 

= {viD — a/\/2) + a/V2, r]{D — al\f2) + ajV2j . 

The bounding curves of the remaining blocks that constitute the 

square with the cavity at the center are obtained via rotation by a factor 7r/2, 


rW rcos7r/2 -sin7r/2 

b *-ii sin7r/2 cos7r/2 


* = 2,3,4, j = S,N,W,E. 


(18) 


The domain can now be represented by attaching the grid representing the 
cavity to one or more Cartesian blocks. 

( 2 ) ( 2 ) 

The circular inclusion of 222 is decomposed into five blocks Bl -B^ . The 
( 2 ) 

block 13]^ Ms a square at the center of the circular inclusion with corners at the 
points (± 0 ( 2 , ±ad), 0 < d < \/2l2 defined by the bounding curves, 


= a {2di - d,-d ), = a {2d^ -d,d), 

= a {—d, 2 ( 2*7 ~ d ), = a {d, 2dr] — d). 

( 2 ) 

Here a is the radius of the circular inclusion. The block B 2 is defined by its 
bounding curves. 


C 

C 

C 


( 2 ) 

2S 


( 2 ) 

2N 


( 2 ) 

2W 


^(2) 

'-'2E 


_ (-( 2 ) 

_ ^(1) 

= a (^—r]{V2/2 — d) — d, rj{-\/2/2 — d) + ctj 
= a (^ri{'/2/2 — d) + d, r\{\f2l2 — (2) ± (2^ . 


( 2 ) ( 2 ) 

The bounding curves of the remaining blocks B\ -B^ that constitute the cir¬ 
cular inclusion of 2?2 are obtained via rotations by a factor 7r/2 as in (TTSl) . The 
inclusion is then centered inside a square of side 2D by attaching it to the blocks 
b[^^-B2^^ above. The circular inclusion is now the union of the nine blocks b[^^- 
B^^ and B^^^-B^\ The domain 2?2 can now be represented by attaching the 
grid representing the inclusion to one or more Cartesian blocks. 
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